added a tiny little helper to prevent zero division and running into errors

Will take a look as soon as I get time.

Couple of notes and questions:

• What are the circumstances that leads to division by zero? Would it be better to mitigate beforehand?
• Why is division by zero a problem here? I assume values become inf. Would it be better to mitigate afterwards? You could easily replace inf values in the output. This is more flexible than clipping output inside eig_special_3d.
• If we clip values to avoid div by zero, we should use probably use something like numpy.nextafter.
• The same solutions needs to be applied in st3dcp.py.

Sorry for my late reply, and thanks very much for pulling together your questions. To your comments:

• In my case the µCT data has relatively low contrast and a voxel size of only 2.3 µm/voxel for carbon fibers (~7 µm diameter, example), meaning the fibers are only partially resolved. As a result, image gradients can be weak and noisy, which makes the gradient estimation locally unstable. In such regions the structure tensor becomes locally isotropic (i.e. eigenvalues become very similar or very small). This can lead to very small or near-zero normalization factors, and consequently to division-by-zero during vector normalization. Mitigating this beforehand avoids propagating inf/nan values into subsequent computations.

• Yes — division by (near) zero produces inf/nan in the normalized vectors, and those can leak into downstream computations (e.g. orientation statistics or cosine similarity). Therefore I think it is preferable to avoid generating them in the first place, rather than cleaning them up afterwards.

• I also agree with your suggestion that if we introduce a numerical floor, using something like numpy.nextafter is preferable to a hard-coded constant. For example:

eps = np.nextafter(
    np.array(0, dtype=l.dtype),
    np.array(1, dtype=l.dtype)
)
l = np.maximum(l, eps)
vec /= l

This guarantees a strictly positive minimum representable value for the corresponding dtype and avoids introducing an arbitrary threshold like 1e-24. Also expierenced vector magnitudes >1 when added a threshold like clipping.

• For the CuPy version, the same approach should apply equivalently using the backend (cp.nextafter). Or something like lib.nextafter depending on the availability and import.

No problem. Sorry for slow reply too. Been on vacation.

While I understand that it's frustrating getting back values like inf, I'm not convinced that it's the job of the eig_special_3d or the structure tensor library to try and "fix" numerical precision issues by replacing infinite (or NaN) values.

My main arguments against are:

1. Concerns I don't think it should be the concern of eig_special_3d to guard against numerical precision issues. It's up to the caller of the function to either provide inputs that result in numerically stable outputs, or clean outputs. I think this is how most numerical libraries like Scipy behave.
2. No correct way There's not one correct way to clean/clip numerical errors. However we chose to do it would only be "correct" to some users.
3. Side effects/unexpected behavior Clipping/replacing specific values, especially in way users don't expect, may lead to side effects that users are not aware of.
4. Cleaning output is simple and flexible The user can easily replace inf values in the way that fits their needs using something like vec[vec == np.inf] = .... Likewise, the caller can easily see which values we effected. If eig_special_3d replaces inf values, the caller has no information about if or where numerical issues actually occurred.

Perhaps I'm missing something here? Is there a reason you can't simply clean the outputs of eig_special_3d?

One more thing. Is the input float32 or float64. If it's 32-bit, I'd recommend trying out 64-bit instead to reduce the risk of numerical errors.

Thanks, that makes sense — and I think I understand your concern more clearly now.

I agree that silently applying a cleanup policy inside eig_special_3d can be problematic, especially if users expect full visibility into where degeneracies occurred and may want to handle them differently downstream.

My main motivation for the pull request is that in structure-tensor applications, near-isotropic regions (e.g., matrix-rich areas in composites) are not invalid input or just floating-point noise. They naturally occur in homogeneous or very low-gradient regions. In these cases the structure tensor itself can approach the zero matrix, since it is computed from local image gradients. As a result, all eigenvalues can become very small and the eigenvector direction becomes ill-defined, so normalization can hit zero or near-zero magnitudes. From my perspective this feels closer to handling an algorithmic degeneracy of the structure tensor rather than purely post-hoc output cleaning.

Currently I am working with float32, but I believe this is not fundamentally a precision issue. Using float64 may reduce the likelihood of hitting exact zeros, but degenerate or near-zero structure tensors will still occur in homogeneous regions — it would mainly shift the threshold at which the problem appears rather than eliminate it.

They naturally occur in homogeneous or very low-gradient regions. In these cases the structure tensor itself can approach the zero matrix, since it is computed from local image gradients. As a result, all eigenvalues can become very small and the eigenvector direction becomes ill-defined, so normalization can hit zero or near-zero magnitudes. From my perspective this feels closer to handling an algorithmic degeneracy of the structure tensor rather than purely post-hoc output cleaning.

Yes, I see that. The question is what behavior we expect of the functions, i.e., what is "correct" behavior and what's considered valid input and output. Right now we have two parts:

1. Calculate the structure tensor, $S$, using structure_tensor_3d.
2. Calculate the eigenvectors and values using eig_special_3d.

The current train of thoughts go something likes this: Any volume of numbers values should be valid input for structure_tensor_3d. So, a volume of zeros is valid input. The output would be, $S$, consisting of just zeros. Then eigenvalues are 0 and eigenvectors inf. This makes sense, as we cannot calculate the eigenvector for a zero-vector.
The same thing may happen to near-zero values of $S$. To me that seems correct. The output indicates there's too little structure to calculate a direction at the given precision. In that case, the user can choose to use higher precision, or just accept that there's practically no structural information and deal with is as they see fit.

We could change inf to 0 or nan, but is there a good reason to do so?

Currently I am working with float32, but I believe this is not fundamentally a precision issue. Using float64 may reduce the likelihood of hitting exact zeros, but degenerate or near-zero structure tensors will still occur in homogeneous regions — it would mainly shift the threshold at which the problem appears rather than eliminate it.

If you actually have structure in all regions of the data, I think float64 will likely solve the issues you've had. This is probably the simplest and most accurate solution.

@choROPeNt as of now, I'm not in favor of clipping values. I really appreciate your input though. I think the best approach is usually to use float64 for calculations to avoid numerical issues. This is also the reason I always cast to float64 in the multiprocessing code.

In any case, it is probably a good idea to very/clean output as needed. I prefer to leave that to the caller, as they know best how to handle numerical instability in their data.